************************************************************************************************
* This replication file estimates a free entry counterfactual with a potential price reduction *
************************************************************************************************
*1. Re-calibrate parameters to reflect price reduction.
*2. Estimate and store information on the status quo for comparison.
*3. Run the counterfactual to estimate the change in welfare after the price reduction 
*4. Store the outcomes  


*************************
* 0 Set path to files *
************************* 
clear all
capture log close 
set more off
set trace off
mata: mata set matastrict off
cd $main_directory


**************************************************************************************
* 1. Load the parameters and re-calibrate the intercept to reflect price reductions. *
**************************************************************************************

use "Data_Belgium\Generated_data\baseline_demand_dataset.dta", clear
qui do "Command_files\spatial_demand_Q.do"

//Set predetermined values:
sca scale=10
set seed 1234
mata: mata matuse "Estimates\betas_Q", replace
mata: mata matuse "Estimates\betas_Q1", replace
mata: mata matuse "Estimates\betas_Q2", replace
mata: mata matuse "Estimates\cost_price_parameters", replace
mata:    st_numscalar("p_Q1",p_Q1)
mata:    st_numscalar("p_Q2",p_Q2)
mata:    st_numscalar("alpha",alpha)
mata:    st_numscalar("nu",nu)
mata:    st_numscalar("aL",aL)
mata:    st_numscalar("phi",phi)
mata:    st_numscalar("aM",aM)
mata:    st_numscalar("fixed_cost",fixed_cost)

//Set predetermined values:
//global reduction_level_Q1 35 //This is the percentage change in the price of real-estate transactions if you would like to set it manually.
//global reduction_level_Q2 26 //This is the percentage change in the price of other transactions if you would like to set it manually.
mata: change_p_Q1=1-`=$reduction_level_Q1'/100
mata: change_p_Q2=1-`=$reduction_level_Q2'/100
global starting_value `" 1 "' //If this is set to 1, then we begin with the status quo starting values.

mata: lambda=1
mata: profit_constraint=1

//Reflect the change in the price:
mata: p_Q1_new=change_p_Q1*p_Q1 //Average price
mata: p_Q2_new=change_p_Q2*p_Q2 //Average price

//Adjust the parameter estimates to take into account the new price
mata: betas_Q1_new=betas_Q1
mata: betas_Q2_new=betas_Q2 
mata:    st_numscalar("p_Q1_new",p_Q1_new)
mata:    st_numscalar("p_Q2_new",p_Q2_new)
//beta_star reflects how much a Euro of transportation costs is worth in terms of utility
mata: beta_star_Q1=betas_Q1[1]/transportation_cost //The transportation costs can be either in thousands or in other units - this does not matter, as long as the price has a similar term of measurement.
mata: beta_star_Q2=betas_Q2[1]/transportation_cost //The transportation costs can be either in thousands or in other units - this does not matter, as long as the price has a similar term of measurement.
//gamma_star reflects the intercept net of price effects
mata: gamma_star_Q1=betas_Q1[rows(betas_Q1)-1]-beta_star_Q1*p_Q1
mata: gamma_star_Q2=betas_Q2[rows(betas_Q2)-1]-beta_star_Q2*p_Q2
//The new intercept is gamma_star plus the new price effects
mata: betas_Q1_new[rows(betas_Q1_new)-1]=gamma_star_Q1+beta_star_Q1*p_Q1_new
mata: betas_Q2_new[rows(betas_Q2_new)-1]=gamma_star_Q2+beta_star_Q2*p_Q2_new

***********************************************************************
* 2. Estimate and store information on the status quo for comparison. *
***********************************************************************
//Calculate the total demand for real estate on the market
qui su Q1 if common==1
sca sumQ1=r(sum)
sca di sumQ1
//Calculate the total demand for other transactions on the market
qui su Q2 if common==1
sca sumQ2=r(sum)
sca di sumQ2
//Calculate the total population on the market
bysort market_ID:  gen dup_market = cond(_N==1,0,_n)
qui su population if dup_market<2
sca sumPopulation=r(sum)
sca di sumPopulation
//Calculate total potential demand per capita as 10 times larger than the current value.
scalar gamma_Q1=scale*sumQ1/sumPopulation
sca di gamma_Q1
scalar gamma_Q2=scale*sumQ2/sumPopulation
sca di gamma_Q2
//Calculate the market size (i.e. number of possible transactions)
gen market_size_Q1=gamma_Q1*population
gen market_size_Q2=gamma_Q2*population

//Define the product for which demand will be estimated.
gen firm_demand_Q1=Q1 
gen firm_demand_Q2=Q2

//Generate a unique numeric identifier for the market.
sort market_ID
egen market_ID_num = group(market_ID)

// Controls
global controls_utility distance lnNotPerOff start_not LDE Perc_Young Perc_Elderly log_income popdens
global controls_total_Q1 newid firm_demand_Q1 market_ID_num market_size_Q1 $controls_utility
global controls_total_Q2 newid firm_demand_Q2 market_ID_num market_size_Q2 $controls_utility

//Define newid, which is just firm_IDs which go from 1 to 1152.
egen newid = group(firm_ID)
//Generate a unique observation identifier, which is used to merge the results from the counterfactuals with the initial dataset
tostring market_ID_num, gen(market_ID_string)
tostring newid, gen(newid_string)
gen join="_"
egen unique_identifier=concat(newid_string join market_ID_string)
drop join market_ID_string newid_string

// Count the number of firms in the sample
bysort firm_ID:  gen dup_firm = cond(_N==1,0,_n)
count if dup_firm<2
sca n_firm=r(N)
// Count the number of markets
bysort market_ID: gen nvals = _n == 1 
count if nvals
sca markets_number=r(N)
di markets_number
drop nvals
* Count the number of markets with a local notary
preserve
keep if common==1 //This denotes the local market of a notary
bysort market_ID: gen nvals = _n == 1 
count if nvals
sca covered_markets_current=r(N)
sca uncovered_markets_current=markets_number-covered_markets_current
su NotPerOff
sca total_notaries_current=r(sum)
restore

//Status quo estimation
* Sort the data by market ID, since this will be the first relevant grouping (to calculate notary utility within a municipality)
sort market_ID_num newid
mata:
// Select the relevant variables and read them into mata
X_Q1    = st_data(., "$controls_total_Q1")
X_Q2    = st_data(., "$controls_total_Q2")
// nx counts the number of observations
nx    = rows(X_Q1)
// Add an intercept, with length equal to the number of observations and 
X_Q1    = X_Q1,J(nx,1,1)
X_Q2    = X_Q2,J(nx,1,1)

// Calculate per firm profit as a matrix 
n_firm=st_numscalar("n_firm")
total_notaries_current=st_numscalar("total_notaries_current")
Results_ln_current_Q1=spdemand_log(X_Q1,betas_Q1)
X_Q1=sort(X_Q1,(3,1)) //sort by market and then firm_id
Results_ln_current_Q2=spdemand_log(X_Q2,betas_Q2)
X_Q2=sort(X_Q2,(3,1)) //sort by market and then firm_id
Results_current_Q1=exp(Results_ln_current_Q1[,2])
Results_current_Q2=exp(Results_ln_current_Q2[,2])
Results_current=Results_current_Q1+Results_current_Q2
output_current_Q1=round(sum(Results_current_Q1))
output_current_Q2=round(sum(Results_current_Q2))
output_current=round(sum(Results_current))
consumer_surplus_m_current_Q1=consumer_surplus(X_Q1,betas_Q1,transportation_cost)
consumer_surplus_m_current_Q2=consumer_surplus(X_Q2,betas_Q2,transportation_cost)
consumer_surplus_current_Q1=round(sum(consumer_surplus_m_current_Q1[,2]))
consumer_surplus_current_Q2=round(sum(consumer_surplus_m_current_Q2[,2]))
consumer_surplus_market_current=consumer_surplus(X_Q1,betas_Q1,transportation_cost)+consumer_surplus(X_Q2,betas_Q2,transportation_cost)
per_firm_variable_profit_current=per_firm_profit_mult_v2(X_Q1,X_Q2,betas_Q1,betas_Q2,p_Q1,p_Q2,nu,aL,phi,aM,alpha)
consumer_surplus_current=round(sum(consumer_surplus_market_current[,2]))
producer_surplus_current=round(sum(per_firm_variable_profit_current[,2])-fixed_cost*total_notaries_current)
welfare_current=round(consumer_surplus_current_Q1+consumer_surplus_current_Q2+producer_surplus_current)
    st_numscalar("output_current", output_current)
   st_numscalar("output_current_Q1", output_current_Q1)
   st_numscalar("output_current_Q2", output_current_Q2)
   st_numscalar("consumer_surplus_current_Q1", consumer_surplus_current_Q1)
   st_numscalar("consumer_surplus_current_Q2", consumer_surplus_current_Q2)
   st_numscalar("producer_surplus_current", producer_surplus_current)
   st_numscalar("welfare_current", welfare_current)
//Generate a vector with the number of notaries in each office for every office-market combination
n_notaries=J(nx,1,0)
n_notaries=round(exp(X_Q1[.,6]))
//Generate a vector with the number of notaries for every office
notary_vector=J(n_firm,1,0) //Make an empty notary vector
//Calculate the actual notaries per office under the status quo
for(f=1; f<=n_firm; f++) {
notary_vector[f,1]=sum(n_notaries:*(X_Q1[.,1]:==f))/(sum((X_Q1[.,1]:==f))) //This is essentially calculating the average number of notaries across all observations for the same notary office, which is just the number of notaries for that office.
}
	end

******************************************************************************************
* 3. Run the counterfactual to estimate the change in welfare after the price reduction. *
******************************************************************************************

* Sort the data by market ID, since this will be the first relevant grouping (to calculate notary utility within a municipality)
sort market_ID_num newid
mata:
// Select the relevant variables and read them into mata
X_Q1    = st_data(., "$controls_total_Q1")
X_Q2    = st_data(., "$controls_total_Q2")
// nx counts the number of observations
nx    = rows(X_Q1)
// Add an intercept, with length equal to the number of observations and 
X_Q1    = X_Q1,J(nx,1,1)
X_Q2    = X_Q2,J(nx,1,1)
//Calculate output
Results_ln_current_Q1L=spdemand_log(X_Q1,betas_Q1_new) //The letter L is used to signify the outcomes after lowering the price.
Results_ln_current_Q2L=spdemand_log(X_Q2,betas_Q2_new)
Results_current_Q1L=exp(Results_ln_current_Q1L[,2])
Results_current_Q2L=exp(Results_ln_current_Q2L[,2])
output_current_Q1L=round(sum(Results_current_Q1L))
output_current_Q2L=round(sum(Results_current_Q2L))
Results_currentL=Results_current_Q1L+Results_current_Q2L
output_currentL=round(sum(Results_currentL))
//Calculate profits per firm as well as industry profits
per_firm_variable_profit_curL=per_firm_profit_mult_v2(X_Q1,X_Q2,betas_Q1_new,betas_Q2_new,p_Q1_new,p_Q2_new,nu,aL,phi,aM,alpha) //"curL" stands for "currentL". There is a limit to the number of characters in the names of mata matrices.
producer_surplus_currentL=round(sum(per_firm_variable_profit_curL[,2])-fixed_cost*total_notaries_current)
//Calculate consumer surplus
consumer_surplus_m_current_Q1L=consumer_surplus(X_Q1,betas_Q1_new,transportation_cost)
consumer_surplus_m_current_Q2L=consumer_surplus(X_Q2,betas_Q2_new,transportation_cost)
consumer_surplus_current_Q1L=round(sum(consumer_surplus_m_current_Q1L[,2]))
consumer_surplus_current_Q2L=round(sum(consumer_surplus_m_current_Q2L[,2]))
welfare_currentL=round(lambda*consumer_surplus_current_Q1L+lambda*consumer_surplus_current_Q2L+producer_surplus_currentL)

//Make a vector with the starting number of notaries on the market. If starting value=1, then we start at the status quo. Otherwise, we start of zero entrants.  
n_notaries=J(nx,1,0)
if ($starting_value == 1) {
n_notaries=round(exp(X_Q1[.,6]))
notary_vector=J(n_firm,1,0) //Make an empty notary vector
//Calculate the actual notaries per office right now
for(f=1; f<=n_firm; f++) {
notary_vector[f,1]=sum(n_notaries:*(X_Q1[.,1]:==f))/(sum((X_Q1[.,1]:==f)))
}
}
firm_profitL=per_firm_variable_profit_curL[.,2]-fixed_cost*notary_vector //This is the vector of current profits, which will be updated with each round of added notaries.
unprofitable_firms=sum((firm_profitL:<0))
st_numscalar("unprofitable_firms", unprofitable_firms)
  st_numscalar("output_currentL", output_currentL)
   st_numscalar("output_current_Q1L", output_current_Q1L)
   st_numscalar("output_current_Q2L", output_current_Q2L)
   st_numscalar("consumer_surplus_current_Q1L", consumer_surplus_current_Q1L)
      st_numscalar("consumer_surplus_current_Q2L", consumer_surplus_current_Q2L)
   st_numscalar("producer_surplus_currentL", producer_surplus_currentL)
   st_numscalar("welfare_currentL", welfare_currentL)
	end

sca change_CS_Q1=consumer_surplus_current_Q1L-consumer_surplus_current_Q1
sca change_CS_Q2=consumer_surplus_current_Q2L-consumer_surplus_current_Q2
sca change_producer_surplus=producer_surplus_currentL-producer_surplus_current
sca change_welfare=welfare_currentL-welfare_current
sca change_Q1_perc=round((output_current_Q1L-output_current_Q1)*100/output_current_Q1,0.01)
sca change_Q2_perc=round((output_current_Q2L-output_current_Q2)*100/output_current_Q2,0.01)
mata: change_Q1_perc=st_numscalar("change_Q1_perc")
mata: change_Q2_perc=st_numscalar("change_Q2_perc")
mata: change_CS_Q1=st_numscalar("change_CS_Q1")
mata: change_CS_Q2=st_numscalar("change_CS_Q2")
mata: change_producer_surplus=st_numscalar("change_producer_surplus")
mata: change_welfare=st_numscalar("change_welfare")

**************************
* 4. Store the outcomes. * 
**************************
mata: results_`=$reduction_level_Q1'_`=$reduction_level_Q2'=`=$reduction_level_Q1',`=$reduction_level_Q2',change_Q1_perc,change_Q2_perc,change_CS_Q1,change_CS_Q2,change_producer_surplus,change_welfare,unprofitable_firms,producer_surplus_currentL
mata: mata matsave "Estimates\Price_reduction\results_`=$reduction_level_Q1'_`=$reduction_level_Q2'" results_`=$reduction_level_Q1'_`=$reduction_level_Q2', replace

//Counterfactual outcome reported in Table 10
di producer_surplus_currentL
sca change_CS_Q1=round(change_CS_Q1/1000)
sca change_CS_Q2=round(change_CS_Q2/1000)
sca change_producer_surplus=round(change_producer_surplus/1000)
sca change_welfare=round(change_welfare/1000)
di $reduction_level_Q1 " & " $reduction_level_Q2 " & " change_Q1_perc " & " change_Q2_perc " & " change_CS_Q1 " & " change_CS_Q2 " & " change_producer_surplus " & " change_welfare " & " unprofitable_firms

